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1.  INTRODUCTION 


The  flowing  medium  in  a  gun  tube  typically  is  a  mixture  of  a  compressible  gas  and  burning 
solid  propellant  grains.  Details  of  the  flow  are  important  for  weapons  development,  but  only 
bulk  properties  can  be  routinely  measured,  such  as  the  trajectory  of  the  projectile,  the  pressure 
history  at  a  fixed  station,  the  heating  inside  the  gun  tube,  etc.  Therefore,  a  need  exists  for  a 
detailed  mathematical  model  of  interior  ballistics  two-phase  flows,  and  an  algorithm  to  solve 
the  correspondingequations. 

The  three-dimensional  mathematical  model  is  developed  carefully  in  Reference  [1].  Refer¬ 
ences  [2]  and  [3]  are  shorter  versions  of  Reference  [1].  This  two-phase  model  is  based  on 
instantaneous,  finite  volume,  weighted  averaging,  and  consists  of  nonlinear  partial  differential 
equations,  constitutive  laws  for  the  averaged  variables,  and  correlations  for  the  interphase 
terms.  The  transient  phenomena  included  in  this  model  are:  the  convection  of  the  phases 
driven  by  gas  phase  pressure,  gas  phase  viscous  stresses,  turbulence,  intergranular  stresses,  inter¬ 
phase  drag,  and  interphase  mass  transfer  due  to  burning  of  the  solid  grains;  the  change  of 
energy  in  the  gas  due  to  convection,  pressure,  laminar  and  turbulence  dissipation,  conduction, 
and  interphase  heat  transfer;  and  the  change  of  the  geometry  and  number  of  the  burning  grains. 
These  phenomena  occur  within  the  volume  defined  by  the  gun  tube  and  the  base  of  an 
accelerating  projectile.  The  motion  of  the  projectile  and  the  two-phase  flow  field  are  coupled 
via  the  gas  pressure  exerted  on  the  projectile’s  base. 

This  model  is  specialized  to  the  case  of  axial  symmetry  of  the  flow  within  the  tube.  Refer¬ 
ence  [1]  lists  all  the  differential  equations,  constitutive  laws  and  correlations.  This  axial  sym¬ 
metric  model  and  numerical  scheme  are  encoded  in  the  DELTA  computer  code.  The  purpose 
of  this  paper  is  to  describe  in  some  detail  the  numerical  algorithm  (Section  2),  and  to  present 
some  validating  computer  runs  of  the  combined  model  and  algorithm  Section  3). 

Previous  multi-dimensional,  multi-phase  work  applied  to  interior  ballistics  includes  that  from 
Paul  Gough  Associates,  Inc.  [4-5],  and  Scientific  Research  Associates,  Inc.  [6-7].  Gough’s  work 
addresses  the  inviscid  flow  during  the  ballistic  cycle,  i.e.  the  average  gas  phase  viscous  stresses, 
heat  conduction,  and  turbulence  are  excluded.  Thus,  the  phenomenology  of  pressure  waves 
inside  a  gun  tube  is  modeled  primarily  in  Gough’s  work.  Scientific  Research  Associates  con¬ 
sidered  the  viscous  flow  phenomena,  i.e.  the  development  of  boundary  layers  inside  a  gun  tube. 
The  work  reported  here  resembles  that  of  Scientific  Research  Associates;  however,  the  model  is 
somewhat  different  (the  governing  equations,  choice  of  dependent  variable,  some  correlations) 
and  numerical  algorithm  has  been  modified.  Differences  and  similarities  in  the  models  are 
addressed  in  detail  in  Reference  [1],  and  in  the  algorithm  in  Section  2  of  this  paper. 
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2.  ALGORITHM 


2.1  Governing  Equations 

The  governing  equations  are  a  set  of  nonlinear  partial  differential  equations  which  are  first 
order  in  time  and  second  order  in  the  two  spatial  coordinates.  The  rationale  for  the  specific 
form  of  the  equations  is  given  in  Reference  [1].  The  general  form  is 

y,  =  G(r,z,t,y,yr,yz,yn,yrz,yu),  (2.1) 


where  the  independent  variables  of  time,  radial  coordinate  and  axial  coordinate  are  denoted  by 
t,  r,z,  respectively.  The  vector  of  dependent  variables  is  denoted  by  y,  and  the  partial  deriva¬ 
tives  of  y  with  respect  to  the  spatial  and  temporal  coordinates  are  denoted  by  subscripts.  The 
components  of  the  vector  y  can  be  the  radial,  circumferential  (swirl),  axial  components  of  the 
gas  phase  velocity  ( u,v  ,m>),  respectively;  the  radial,  circumferential  (swirl),  axial  components  of 
the  solid  phase  velocity  (u*  ,v*  ,w’),  respectively;  the  gas  phase  specific  entropy  .f  ;  the  loga¬ 
rithm  of  gas  phase  pressure <7  ;  the  regression  distance  of  the  solid  phase  d'  ;  the  number  of 
particles  m  *  ;  the  particle  surface  temperature  T*  ;  and  two  variables  to  define  the  turbulence 
in  the  flow  field.  Thus,  depending  on  the  simulation,  the  number  of  dependent  variables  change 
from  a  minimum  of  four  to  a  maximum  of  thirteen.  For  the  case  of  an  one-phase,  laminar 
flow  simulation  with  swirl,  the  dependent  vector  y  has  five  components,  u ,  v,  w,  s  and  q.  For 
the  case  of  a  two-phase,  turbulent  flow  simulation  with  ignition  and  burning  of  the  solid  phase, 
the  dependent  vector  y  has  nine  components,  u,  w,u *,  w' ,  s ,  q ,  d* ,  m' ,  T*.  This  assumes  the 
absence  of  swirl  and  an  algebraic  turbulence  model  (i.e.  the  turbulence  properties  are 
described  only  by  algebraic  relations).  The  components  of  the  vector  G  are  nonlinear  functions 
which  can  depend  on  the  variables r,  z,  t,  y,  yr,  yz,  yn,  yn,  yu. 

The  spatial  domain  is  a  confined  volume  within  a  tube  bounded  in  length  by  a  stationary 
wall  (the  gun  breech)  and  the  base  of  an  accelerating  projectile.  The  radius  of  the  tube  can 
depend  on  the  axial  position  from  the  breech  which  is  denoted  by  zB .  The  radial  coordinate 
varies  from  the  axis  of  symmetry  to  the  tube  wall.  The  radial  positions  of  the  axis  and  wall  are 
denoted  by  rA  and  rw  (2),  respectively.  The  axial  coordinate  varies  from  the  breech  to  the  base 
of  the  projectile.  The  axial  position  of  the  breech  and  projectile  base  may  have  a  radial  depen¬ 
dence  which  is  denoted  by  zB  ( r )  and  zp  ( r ),  respectively. 

The  projectile  is  assumed  to  move  as  a  rigid  body.  The  unsteady  projectile  motion  is 
governed  by  the  following  equations: 
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where  wp,  vp,  mp  denote  the  axial  velocity,  circumferential  velocity  and  mass  of  the  projectile, 
respectively.  The  forces  that  retard  the  motion  of  the  projectile  are  those  due  to  air  resistance, 
friction  between  the  projectile  and  tube  wall,  and  gas  leakage  around  the  projectile,  and  are 
denoted  by  Fj ,  FD,  and  Fp,  respectively.  These  retarding  forces  are  assumed  to  be  known  func¬ 
tions.  If  the  tube  is  rifled,  an  additional  phenomenon  is  present  which  causes  the  projectile  to 
rotate,  and  its  mass  to  be  effectively  increased  via  equation  (2.4).  In  this  case  the  angle  of 
rifling  0R  is  nonzero,  and  the  moment  of  inertia  of  the  projectile  Im  must  be  given.  Because 
the  pressure/?  is  determined  from  the  solution  of  governing  equations  of  the  flow  field,  which 
depends  on  the  value  of  wp,  equations  (2.1)  -  (2.5)  represent  a  coupled  system  with  a  moving 
boundary. 


2.2  Numerical  Algorithm 

We  want  to  compute  by  finite  difference  approximations  the  transient  values  of  the  vari¬ 
ables  which  describe  the  fluid  dynamics  of  the  flow  in  the  region  confined  by  the  inner  tube 
wall,  breech,  and  moving  projectile.  One  way  to  calculate  in  this  expanding  computational 
region  is  by  an  “accordion”  type  grid  in  the  axial  direction,  i.e.  the  first  and  last  axial  grid 
points  are  attached  to  the  breech  and  projectile,  respectively,  and  the  mesh  expands  as  the  pro¬ 
jectile  accelerates  down  the  tube.  Thus,  the  physical  grid  moves  in  accordance  with  the  projec¬ 
tile  motion.  With  regard  to  the  spatial  finite  difference  approximation,  the  goal  is  to  obtain  an 
accurate  approximation  to  the  actual  physical  happening.  It  can  be  shown  that  the  finite  differ¬ 
ence  approximations  to  the  physical  variables  in  the  physical  mesh  are  the  same  whether  one 
directly  differences  on  the  physical  grid,  or  one  differences  on  a  transformed  grid  and  then 
transforms  back  to  the  physical  grid.  Higher  accuracy  in  a  transformed  space  is  not  meaningful 
if  it  is  lost  in  the  transformation  back  to  the  physical  space.  Furthermore,  our  physical  grids 
will  be  orthogonal  or  nearly  orthogonal.  Thus,  we  choose  to  compute  finite  difference  approxi¬ 
mations  to  spatial  derivatives  on  the  physical  grid.  An  additional  advantage  of  our  method  is 
that  the  governing  equations  need  not  be  transformed,  and  thus  are  simpler  to  understand  and 
change  in  the  computer  code. 
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The  finite  differencing  of  the  time  derivatives  can  be  of  two  generic  types:  implicit  or  expli¬ 
cit  (See  Reference  [8]).  The  principal  advantage  of  an  implicit  scheme  is  its  superior  stability 
properties  compared  to  an  explicit  scheme.  For  convection-diffusion  type  problems  like  those 
given  by  equation  (2.1),  an  explicit  finite  difference  method  has  two  stability  conditions,  the 
Courant-Friedrichs-Lewy  condition  and  the  viscous  stability  limits.  In  one  dimension,  these 
condition  are: 

At  =  CFL  — -  ,  (2.6) 

c+w  v  ' 

At  =  VSL^-(Ax)2  , 

2 

where  Ax  is  the  spatial  mesh  increment,  c  is  the  sound  speed,  w  is  the  gas  velocity,  p  is  the 
viscosity  and  p  is  the  density.  The  constants  CFL  and  VSL  are  less  than  or  equal  to  one.  For 
simulations  which  involve  boundary  layers,  small  grid  sizes  are  necessary.  Thus,  for  this  type  of 
simulation,  the  time  step  (the  size  of  At)  must  be  proportional  to  the  square  of  the  smallest 
grid  increment  for  an  explicit  scheme.  On  the  other  hand,  most  implicit  schemes  have  no 
corresponding  stability  conditions,  and  significantly  larger  time  steps  based  on  accuracy  con¬ 
siderations  rather  than  stability  can  be  used.  The  basic  disadvantage  of  implicit  algorithms  is 
that  they  tend  to  be  more  complicated  than  explicit  schemes,  and  thus  more  difficult  to  under¬ 
stand  and  implement.  In  particular,  applying  a  standard  implicit  scheme  to  system  of  equations 
(2.1),  we  obtain  a  nonlinear  system  of  algebraic  equations  in  the  variables  at  the  new  time  level. 
This  system  can  be  quite  large  and  complex  because  it  possesses  an  equation  for  each  depen¬ 
dent  variable  and  for  each  grid  point  in  the  two-dimensional  computational  mesh.  These  equa¬ 
tions  are  coupled  via  the  spatial  derivatives.  Iterative  methods  are  the  most  common  solution 
procedure,  but  they  can  be  quite  complex  and  time-consuming  for  such  a  general  system.  To 
mitigate  these  undesirable  characteristics,  we  apply  a  method  developed  in  References  [9-11]. 
A  salient  feature  of  this  method  is  the  temporal  linearization  of  the  nonlinear  terms  to  within 
the  local  truncation  error  of  the  finite  difference  approximation  of  the  time  derivative.  The 
resulting  system  can  then  be  represented  in  a  matrix  equation, 

Ay"  =  b  (2.7) 

where  y"  is  the  vector  of  unknown  dependent  variables  at  the  new  time  level,  and  A  and  b  are 
the  matrix  and  vector  of  values  at  the  known  time  level,  respectively.  Furthermore,  the  matrix 
A  can  be  structured  if  we  decompose  the  time-differenced  equation  of  (2. 1 )  into  two  systems  of 
equations,  each  of  which  involves  the  spatial  derivatives  of  the  unknown  variable  in  only  one 
coordinate  direction.  This  decomposition  or  splitting  is  done  so  that  the  error  incurred  is  of  the 
order  of  the  local  temporal  truncation  error,  and  so  that  the  decomposed  or  split  equations  still 
form  a  consistent  approximation  to  equation  (2.1).  When  centered  differences  are  used  to 
approximate  the  spatial  derivatives,  tridiagonal  matrices  are  obtained.  Because  we  are  dealing 
with  a  system  of  equations,  the  matrices  are  block  tridiagonal  where  the  size  of  the  blocks  is 


equal  to  the  number  of  dependent  variables.  This  method  of  splitting  the  implicitly  differenced 
equations  is  called  an  Alternating  Direction  Implicit  (ADI)  scheme  (see  Reference  [8]). 
Because  we  have  also  linearized  the  equations,  we  shall  refer  to  this  method  as  a  linearized 
ADI  scheme. 

This  general  linearized  ADI  scheme  is  applied  to  the  instantaneous,  finite-volume,  weighted, 
averaged  equations  of  interior  ballistics  with  several  unique  features:  The  algorithm  is  derived 
for  a  moving  coordinate  system,  and  has  no  mass  source  due  to  the  motion  of  the  grid.  The  ele¬ 
ments  of  the  matrices  derived  by  the  linearization  process  are  obtained  directly  by  numerical 
differentiation  which  bypasses  the  tedious  and  error  prone  task  of  analytically  deriving  each  ele¬ 
ment,  and  the  subsequent  coding  of  these  complex  expressions.  The  spatial  differencing  is  per¬ 
formed  directly  on  nonuniform  distributed  grids. 


2.2.1  Algorithm  for  Non- Boundary,  Non-Center-Line  Points 

We  derive  the  numerical  scheme  for  the  system  of  equations  (2.1)  on  a  moving  coordi¬ 
nate  system;  i.e.,  the  coordinates  of  the  spatial  grid  system  varies  in  time.  We  let  the  super¬ 
scripts  n  denote  the  new  time  level  and  c  the  current  time  level.  The  change  in  the  j‘h  coordi¬ 
nate  position  of  a  spatial  coordinate*  from  level  c  to  n  is  denoted  by 

AXj  =  *;-*;  =  (i f"  -tc)  =  O(Ar),  tc  <t*<tn,  (2.8) 

at  t=t' 

where  A t  =  tn  -  tc  .  A  Taylor  expansion  of  y”  =  y (r" ,zn,t")  about  the  current  values  at 
(rc ,  zc ,  tc )  is 

J"  =  f  +  +  Si ±Az  (2.9) 
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By  adding  (1-/3)  times  (2.9)  and  p  times  a  similar  expansion  to  (2.9)  of  yc  expanded  about 
y",  and  noting  that  y"  -  yc  =  O(Af).  z"  -  zc  =  O(At),  and  rn  -  rc  =  O (At),  we  obtain 
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where  ET(At 3)  denotes  the  neglected  truncation  error  of  0(At3).  For  a  stationary  grid 
Ar  =  Az  =  0,  we  obtain  the  standard  Crank-Nicolson  scheme  for  integration  parameter  0  —  2- 

and  the  standard  fully  implicit  scheme  for  0  =  1.0.  Equation  (2.10)  is  a  nonlinear  system  of 
equations  in  y"  because  G"  is  a  nonlinear  function.  To  make  (2.10)  a  system  of  linear  equa¬ 
tions,  we  linearize  G"  via  a  Taylor  expansion  about  the  current  level,  i.e. 


G"  =  Gc  + 
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dt 

k  4 

dr 

J 

dz 

k  > 

Az 


+  DC  Ay  +  DRC  Ayr  +  DZC  Ay, 


+  DRRC  Ayn  +  DRZC  Ayn  +  DZZC  Ay u 


+  Eu(At2), 

where  Eu  (At2)  denotes  the  neglected  linearization  error  of  0(At2).  The  Jacobian  type 
matrices  are  denoted  by  D,  DR,  DZ,  DRR,  DRZ,  DZZ,  are  evaluated  at  the  current  time 
level,  and  are  defined  as 


\dGt  1 

c 

,  DRfi  = 

f9Gi  1 

C 

,  •  •  •  ,  DZZf:  = 

dG, 

dyj 

’  'J 

1  ‘  J 

(2.12) 


Upon  substituting  (2.11)  into  (2.10),  we  obtain  the  following  linear  system  of  equations  in  y" 
after  algebraic  manipulation: 
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y"  -  0  <  Ar  4^  +  Az  ■§■  +  At  Dcyn  +  DRC  +  DZC 

or  dz  or 


+  DRRC 


-^■1  +DRZe 
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dr 2 

drdz 

dz 2 

(2.13) 


yc  -  /3  <Ar  +  A *  41  +  At  Deyc+DRc  41  +  DZC  &L 

i  or  oz  or  dz 


+  DRRC 


+  DRZC 


Jh- 

drdz 


+  DZZC 


ii]C 

dz 2 


'  'ic  >c  f  lc  'ic 

+  A r  +  Az  4**  +A<Gc+j9Af  Ar  4r"  +  Ar  4r- 

or  oz  dt  or 


+  Az  ^  +  (/l-i)0(At2)  +  £r(Af3)  +  £u(Af3). 


The  neglected  linearization  error  Eu  is  now  0(Af3)  because  G"  was  multiplied  by  At  in 
(2.10).  Thus,  the  linearization  process  does  not  alter  the  order  of  the  temporal  error. 

The  values  yc,  rc,zc,  tc ,  rn,zn,  tn  are  all  known  before  the  start  of  the  integration  routine 
to  determine  the  values  of  yn .  Thus,  the  right  hand  side  of  (2.13)  is  a  known  vector.  The  left 
hand  side  of  (2.13)  can  be  written  as  a  matrix  with  known  values  times  a  vector  of  unknowns, 
y" .  Thus,  (2.13)  has  the  form  of  (2.7).  The  matrix  A  interconnects  values  of  yn  at  a  grid  point 
to  all  the  values  of  its  spatial  neighbors  via  the  first  and  second  spatial  partial  derivatives  at 


level  n.  If  the  term  I  4  »  1  were  absent  from  the  left  hand  side  of  (2.13),  we  could  decom- 
drdz  v  ' 


pose  (2.13)  into  two  matrix  equations  which  are  highly  structured  and  easily  solvable  while  still 
retaining  only  two  time  levels  of  the  dependent  variable  vector,  i.e.  y”  and  yc .  To  this  end,  we 
linearize  the  mixed  derivative  at  the  n  level  about  the  current  level  and  retain  only  the  leading 
term: 


’-V*  •>  *  w  * 


The  symbol  I  represents  the  identity  matrix,  and  the  superscript  k  can  be  either  «  ore. 
The  “lagging”  of  the  mixed  derivative  (2.14)  increased  the  error  of  the  approximation  to 
0(At2)  for  any  value  of  the  integration  parameter  0,  but  we  gain  a  structured  matrix. 

We  can  decompose  (2. 15a)  along  coordinate  directions  in  the  following  manner: 


I  -  0D?  y1  =  I  -  0Dcr  yc  +  Lyc  , 

(2.16) 

I  -  0D"]yf  =  [i-  PDczy  +  (y/  -  yc). 

(2.17) 

Equation  (2.16)  constitutes  the  radial  sweep  of  this  Alternating  Direction  Implicit  method 
because  it  involves  only  spatial  derivatives  in  the  radial  direction  at  the  new  time  level.  One 
solves  this  equation  for  each  fixed  axial  index  and  for  radial  indices  varying  from  the  axis  of 


fea 


ft 


a 


i 


symmetry  to  the  gun  tube  wall.  When  three  point  centered  spatial  finite  differences  are  used  to 
approximate  the  spatial  derivatives  in  the  radial  direction,  the  matrix  I  -  0D"  is  a  block  tridiag¬ 
onal  matrix.  The  size  of  each  block  is  equal  to  the  number  of  dependent  variables.  The 
number  of  block  rows  is  equal  to  the  number  of  grid  points  from  the  axis  of  symmetry  to  the 
gun  tube  wall,  denoted  by  JRMX.  The  block  rows  from  the  second  grid  points  to  one  away 
from  the  wall,  namely  JRMX-l  is  determined  by  (2.16).  The  entries  of  the  first  and  last  block 
rows  are  determined  from  the  conditions  imposed  at  the  axis  of  symmetry  and  wall,  respec¬ 
tively.  The  right  hand  side  of  (2.16)  is  a  known  vector  because  it  is  evaluated  at  the  current 
time-step.  The  solution  of  this  equations  is  the  intermediate  values  of  the  dependent  variables 


Equation  (2.17)  constitutes  the  axial  sweep  of  this  two  sweep  scheme  because  it  involves 
only  the  spatial  derivatives  in  the  axial  direction  at  the  new  time  level.  One  solves  this  equation 
for  each  fixed  radial  index  and  for  axial  indices  varying  from  the  breech  to  the  base  of  the  pro¬ 
jectile.  The  matrix  I  -  0DZ  is  a  block  tridiagonal  matrix  when  three  point  centered  spatial  fin¬ 
ite  differences  are  used  to  approximate  the  partial  derivatives  in  the  axial  direction.  The  size  of 
the  blocks  are  the  same  as  in  the  radial  sweep,  and  the  number  of  block  rows  is  equal  to  the 
number  grids  points  placed  from  the  breech  to  the  projectile  base,  denoted  by  JZMX.  Equation 
(2.17)  is  used  to  determine  the  entries  in  the  block  rows  from  the  second  point  (one  after  the 
breech)  to  JZMX -l  (one  before  the  projectile  base).  The  entries  in  the  first  and  last  block 
rows  are  determined  from  the  boundary  conditions  imposed  at  the  breech  and  projectile  base, 
respectively.  Because  the  intermediate  value  y1  is  known,  the  right  hand  side  of  (2.17)  is 
known.  The  solution  of  (2.17)  is  the  value  of  the  dependent  variables  at  the  new  time  level, 
denoted  by  yF .  The  values  of  the  solution  vector  yF  of  (2.17)  and  those  of  the  solution  vector 
y"  of  (2.15)  differ  by  the  time  error  introduced  by  the  splitting  (2.16)-(2.17).  To  determine  the 
order  of  this  error,  we  substitute  (2.17)  into  (2.16)  and  obtain 


(I  -  0Dr)  [(I  -  0DZ) yF  -  (I  -  0Dz)y '  +  yc  ]  =  [1  -  0Dr]yc  +  Ly' 


(2.18) 


which  simplies  to 

(I  -  0Dr  -  0DZ) yF  =  (I  -  0Dr  -  0DZ) yc  +  Lyc  -  0>DrDz{ yF  -  yc\  (2.19) 


Subtracting  (2.15a)  from  (2. 19),  have 

(I  -  0Dr  -  0Dz)tf  -  y")  =  ~02DrDz(yF  -  y”  +  y"  -  y‘) 


(2.20) 


(I  -  0Dr  -  0DZ  +  fiDrDz){ yF  -  y")  =  -^DrDz( y"  -  yc> 


(2.21) 


We  note  that  the  coefficient  of  yF  -  yn  is  0(1),  Dr  and  Dz  are  each  O(At)  and  (y"  -  yc)  is 


JVW  .Y.V  .>  VX  JYVXV  J.  ''  V 


VW'/  '.-  VV’.>  '  - 


•  V  V  V.V  V.V.V r. 


at  least  O (At).  Thus, 


yf  =  y"  +  £$(A'3). 


that  is,  yF  is  equal  to  y"  to  within  the  local  truncation  error  of  the  scheme. 


(2.22) 


Equations  (2.16),  (2.17)  with  the  definitions  (2. 15b)-(2.15d)  represent  the  time  differenced, 
linearized,  ADI  scheme.  We  now  turned  to  the  finite  difference  approximation  of  the  spatial 
derivatives.  The  standard  centered  finite  difference  approximations  to  the  spatial  partial  deriva¬ 
tives  in  the  coordinate  direction*  at  the  irt  grid  point  are 


iil  n  hi  (y‘  ~  y-i)  fy-ifa+i  -  y.) 

.  &  ),•  hi-\{K  +  K+\)  hi  (hi  +  hi+l)  ’ 


(2.23) 


i 

i 


where  hx 


=  Xi+l-Xi. 


dh  „  2(y/-i  ~  y.)  2(yl>1-y,) 

chc2  i  hi_x(hi  +  hi_x)  ^  +  hi_x)  ’ 


However,  instead  of  (2.23)  we  use 

fo,  _  (y«+i  -  y«-i) 

^  (hi  +  hi_ ,)  ■ 


(2.24) 


(2.25) 


(See  Reference  12).  For  equally  spaced  meshes  (2.23)  and  (2.25)  are  identical.  For  a  nonuni¬ 
form  spaced  grid,  the  difference  between  them  can  be  expressed  as 

T(*,-i  -  (2-26) 


The  nonuniform  grids  that  are  used  in  our  application  have  the  property  that  hj_x  >  ht  so  that 
(2.26)  acts  as  a  stabilizing  viscous  term  to  the  spatial  differencing. 

To  complete  the  description  of  the  interior  point  algorithm  we  discuss  the  determination  of 
the  Jacobian  type  matrices  D,  DR,  DZ,  DRR,  DZZ  defined  by  (2.12).  The  obvious  way  to 
evaluate  the  elements  of  these  matrices  is  to  manually  take  the  partial  derivatives  and  code 
each  element.  This  double  procedure  is  very  error  prone  because  the  right  hand  sides  of  the 
equations  are  extremely  complex.  If  there  are  NEQ  variables,  then  each  element  of  the  vector 
G  is  a  function  of  6 *NEQ  +  3  arguments,  in  general.  Furthermore,  for  two-phase  simulations, 
the  correlations  are  not  fixed,  but  can  vary  substantially  from  simulation  to  simulation.  In  these 
cases  new  elements  of  the  matrices  would  have  to  be  determined  and  encoded.  To  avoid  this, 
one  can  lag  the  contribution  of  these  correlations  by  one  time  step.  This  is  not  totally  desirable 
because  one  increases  the  local  truncation  error  (see  the  discussion  near  (2.14))  which  can  be 
large  when  the  correlations  add  significantly  to  the  flow  dynamics  in  a  given  time  step,  e.g., 
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burning  of  the  grains.  An  alternative  to  this  whole  procedure  is  to  determine  these  matrices 

QQ 

numerically.  Consider  the  determination  of  /),.  =  — -  which  can  be  approximated  by 

dy ) 

Dij  «  ,yj  +*,-■■  ,yNEQ,yr,  -  ,Ja)  (2.27) 

—  (r>  Z,  t,  y  i  ,  •  •  ,  )>j  —  6  , '  , y^jEQ  ’  y r*  i  y zr)  / (25)  , 

where  6  is  the  pre-determined  increment.  Once  these  increments  are  obtained,  the  elements  of 
the  matrices  can  be  computed  trivially  by  repeated  calls  to  a  subroutine  which  computes  the 
right  hand  sides  of  the  equations.  A  characteristic  of  the  G 's  is  that  the  terms  yr  and  y2 
appear  at  most  quadratically,  and  the  terms  yIT  and  yB  appear  at  most  linearly.  Thus,  by  using 
centered  differences  the  matrices  DR  and  DZ  can  be  obtained  exactly  for  any  value  of  the 
increment.  Likewise,  using  one-side  differences  the  matrices  DRR  and  DZZ  can  be  deter¬ 
mined  exactly.  In  these  cases,  a  large  value  of  the  increment  can  be  used  to  avoid  any  round 
off  errors.  On  the  other  hand  the  vector  G  is  a  non-algebraic,  nonlinear  function  of  the  vector 
y.  In  this  case  one  cannot  obtain  an  exact  value  of  the  elements  of  D  for  any  value  of  the 
increment.  We  have  developed  a  strategy  to  compute  an  increment  value  based  on  the  current 
error  estimates  of  G,  and  y y,  that  is,  to  determine  a  &.•  is  used.  The  drawback  of  this 
approach  is  the  computing  time  necessary  to  evaluate  the  right  hand  sides  as  often  as  required. 

Finally,  we  address  the  problem  of  artificial  mass  sources  induced  solely  by  the  motion  of  a 
grid  system.  The  problem  was  illustrated  and  resolved  in  Reference  13.  There  exists  a  standard 
procedure  to  determine  if  mass  sources  occur  in  a  numerical  scheme  when  the  grids  are  moved. 
First  one  assumes  a  constant  flow  field  at  the  current  level,  secondly  one  applies  the  method  to 
compute  the  new  time  level  of  values  on  a  displaced  grid,  and  finally  one  determines  if  these 
new  values  differ  from  the  constant  values.  Following  this  procedure,  we  assume  that  the  flow 
field  variables  are  constants  which  satisfy  the  partial  differential  equations,  and  assume  that 
Ar  *  0  and  Az  *  0.  Consequently,  all  the  spatial  derivatives  at  the  current  time  level  are  zero 
and  (2.13) reduces  to 

y"  —  |at  y"  +  Az  y"  (2.28) 

+  A  t  D  yn  +  DR  y"  +  DZ  y"  +  DRR  y*  +  DRZ  y"  +  DZZ  yB 

j 

=  yc  -  fi  At  D  yc. 


11 


Using  the  fact  that  the  spatial  derivativesat  the  current  level  are  zero,  we  add  zero  to  (2.28)  in 
a  convenient  form  to  obtain 


A(y"  -  yc)  =  0, 


(2.29) 


where  A  is  a  matrix.  If  A  is  nonsingular,  then  the  solution  of  (2.29)  is  y"  =  yc .  Thus,  no  mass 
sources  exist.  The  lack  of  mass  sources  is  due  to  the  form  of  the  equations,  i.e.  y.  =  G,  and  the 
algorithm.  For  a  simpler  set  of  equations  than  the  one  we  are  solving,  and  for  a  set  in  a  con¬ 
servation  form  which  are  transformed  to  a  stationary  uniform  computational  grid,  a  “Geometric 
Conservation  Law”  is  needed  to  prevent  mass  sources.  (See  Reference  13).  Our  method 
automatically  avoids  this  other  partial  differential  equation,  and  the  need  to  obtain  its  solution 
at  every  time  step. 


2.2.2  Algorithm  for  Points  at  the  Center- Line  and  at  the  Solid  Surfaces 

The  method  to  obtain  the  new  time  level  of  the  flow  variables  along  the  center-line  is  dif¬ 
ferent.  To  maintain  the  axial  symmetry  of  the  flow,  the  physical  conditions  on  the  flow  are  the 
radial  and  circumferential  velocity  components  for  both  the  gas  and  particle  phases  must  be 
zero,  and  the  first  partial  derivatives  with  respect  to  the  radial  direction  at  the  axis  of  symmetry 
of  the  remaining  variables  must  be  zero. 

The  contribution  from  the  points  on  the  axis  of  symmetry  can  be  done  in  at  least  two  ways. 
The  first  and  simplest  is  to  directly  apply  the  symmetry  conditions.  For  a  radial  sweep,  the  ele¬ 
ments  of  the  first  block  row  of  the  matrix  A  and  known  vector  b  are  the  finite  difference 
approximations  of  these  conditions.  The  axial  sweep  along  the  center-line  is  performed  after 
the  axial  sweeps  along  interior  axial  indices.  The  final  values  at  the  center-line  are  obtained 
using  the  symmetry  conditions,  and  the  final  axial  sweep  values  of  the  non-center-line  points. 
The  second  way  is  more  complex.  Because  the  center-line  is  part  of  the  flow  field  (physically  a 
non-boundary),  the  governing  partial  differential  equations  are  valid  on  the  axis  of  symmetry. 
One  may  rewrite  these  equations  with  the  symmetry  conditions  imposed  in  the  equations  them¬ 
selves,  and  with  the  correct  limit  conditions  as  the  radial  coordinate  goes  to  zero.  Then,  solve 
these  new  equations  by  exactly  the  some  method  as  described  for  non-boundary  points. 
Although  both  are  coded,  the  simpler  first  option  is  utilized. 

The  boundary  conditions  at  the  solid  surfaces  such  as  the  breech,  tube  wall  and  projectile 
base  can  vary  substantially  with  the  particular  simulation.  This  makes  a  general  discussion  of 
boundary  conditions  difficult.  However,  we  will  discuss  some  implementations  of  common 
types  of  boundary  conditions.  The  simple  functional  form  boundary  condition  y  =  constant , 

and  the  simple  derivative  boundary  condition  =  constant  are  the  most  trivial  kinds.  Their 

OX 

finite  difference  approximations  are  straightforward,  if  y  is  a  variable  computed  directly  by  a 
governing  partial  differential  equation,  i.e.,  one  of  the  variables  in  (2.1).  However,  if  one  has  a 
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nonlinear  boundary  condition  of  the  form  /  |y  J "  =  0  or  /  =0,  then  one  must  linear¬ 

ize  function  /  in  time  in  the  same  manner  as  is  any  component  of  the  nonlinear  vector  function 
G  in  (2.1)  and  (2.11). 

Sometimes  a  nonlinear  boundary  condition  can  be  reformulated  as  a  linear  one.  Consider 


the  adiabatic  condition  —  =  0  were  T  is  the  temperature  and  n  is  the  outward  normal. 

In  our  method  the  entropy  s  and  pressure  function  q  are  computed  directly  from  the  governing 
partial  differential  equations.  Thus  T  is  a  nonlinear  function,  T  —  T(s,q).  We  ordinarily 


would  expand  in  a  Taylor  series  in  time  to  obtain  a  linear  function  in  sn  and  qn . 

However,  for  a  Noble-Abel  equation  of  state,  we  have  a  simplification.  We  use  the  chain  rule 
to  obtain 


dT  ds  |  dT  dq  _  Q 

ds  dn  dq  dn 


(2.30) 


Since  ~~  =*  0  ,  then  (2.30)  can  be  written  as 

OS 


Jl  +  i2_  =0L 

dn  dT  dn 
ds 


(2.31) 


By  noting  that  the  ratio  of  partial  derivatives  of  temperature  is  a  constant,  (2.31)  is  a  linear 
function  in  sn  and  qn  Thus,  (2.31)  can  be  finite  differenced  and  incorporated  directly  into  the 
matrices  of  the  linearized  ADI  method. 


2.2.3  The  Order  of  the  Sweeps 

The  order  of  the  sweeps  is  mainly  a  bookkeeping  problem,  and  should  not  have  a  large 
effect  on  the  solution.  We  have  used  the  following  procedure.  First,  radial  sweeps  from  the 
center-line  to  the  tube  wall  are  performed  along  constant  axial  indices  to  obtain  intermediate 
values  of  the  variables  y  using  (2.16).  The  axial  index  varies  from  the  one  after  the  breech  to 
the  one  before  the  projectile  base.  Second,  axial  sweeps  from  the  breech  to  the  projectile  base 
are  performed  along  constant  radial  indices  to  obtain  values  of  the  variables  y  at  the  new  time 
level  using  (2.17).  The  radial  index  varies  from  the  one  after  the  axis  of  symmetry  to  the  one 
before  the  gun  tube  wall.  Third,  the  symmetry  conditions  are  applied  to  determine  the  new 
time  level  values  at  all  points  on  the  axis  of  symmetry  by  using  the  results  of  the  axial  sweep 


are  determined  by  imposing  the  wall  boundary  condition  using,  if  necessary,  the  results  of  the 
axial  sweep  (step  two). 


We  note  that  the  radial  sweeps  along  the  breech  and  projectile  are  avoided.  The  justifica¬ 
tion  is  that  in  many  applications  the  boundary  conditions  at  the  breech  and  projectile  do  not 
involve  radial  derivatives.  Consequently,  the  linearized  version  of  these  conditions  can  be 
expressed  without  the  need  for  the  determination  of  the  intermediate  values  of  the  variables. 
Recall  that  the  matrix  D  (equation  (2.12))  was  included  in  the  radial  sweep.  (See  equations 
(2.15b)  &  (2.16)).  However,  it  could  have  been  incorporated  into  the  axial  sweep  formula. 
(See  equations  (2.15c)  &  (2.17)).  If  no  radial  derivatives  exist  for  the  boundary  conditions,  we 
avoid  the  radial  sweep  and  incorporate  the  D  matrix  in  the  axial  sweep  formulation.  An  exam¬ 
ple  of  this  type  of  boundary  condition  is  the  gas  continuity  equation  used  to  determine  the 
transformed  pressure  <7  at  the  breech  for  a  one-phase,  viscous  simulation.  Imposing  the  no-slip 
conditions  u  =  v  =  w  =  0  in  the  continuity  equation  evaluated  at  the  breech,  all  radial  deriva¬ 
tives  vanish,  and  the  above  method  is  applied  directly  during  the  axial  sweep  where  one-sided 
finite  differences  are  used  to  approximate  the  spatial  derivatives  at  the  breech.  If  one  uses  the 
normal  velocity  equation  for  the  determination  of  q  at  the  projectile  base,  for  example,  some 
radial  derivatives  remain  in  the  viscous  terms.  Only  by  lagging  them  could  one  use  the  above 
approach. 
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3.  RESULTS 


The  results  of  two  simulations  using  the  DELTA  code  are  presented  in  this  section.  These 
particular  calculations  are  selected  because  they  can  be  compared  to  independently  determined 
answers,  and  thus,  give  some  verification  of  the  code’s  accuracy  and  capabilities.  Both  cases 
involve  a  one-phase  gas  expansion  in  a  constant  cross-section  tube  closed  at  one  end  by  a  sta¬ 
tionary  surface  called  the  breech,  and  at  the  other  end  by  a  movable  piston  called  the  projectile. 
The  breech  and  projectile  base  are  assumed  to  be  flat  surfaces.  The  initial  states  of  the  gas  are 
uniform  and  quiescent,  but  the  gas  pressures  are  great  enough  to  accelerate  the  projectile 
through  the  tube.  The  controlling  mechanisms  of  the  expansion  flows  are  the  same,  namely  the 
propagation  of  the  rarefaction  wave,  generated  by  the  projectile  displacement,  and  its  reflection 
from  the  breech,  then  the  projectile,  and  so  forth.  However,  the  gas  pressure  levels  differ 
greatly  between  the  simulations,  and  the  subsequent  flows  are  in  different  regimes. 


TABLE  1.  Geometry  and  Gas  Properties  for  the  150-mm  Gun 
and  Bicen-Whitelaw  Experiments 


150-mm  Gun  Description 

Bicen-Whitelaw 

150.0 

Bore  Diameter  (mm ) 

76.7 

1.698 

Initial  Projectile  Displacement  (m ) 

0.1773 

6.0 

Maximum  Travel  of  Projectile  (m  ) 

0.3 

50.0 

Projectile  Mass  (kg) 

2.54 

0.001 

Covolume  (m3/kg) 

0.0 

1.22 

Ratio  of  specific  heats,  7 

1.4 

621.09 

Initial  Pressure  (MPa ) 

0.28 

2666.8 

Initial  Temperature  (A) 

293.0 

The  first  simulation  corresponds  to  the  gas  expansion  within  a  150-mm  tube  away  from  the 
effects  of  the  tube  wall  (core-flow)  under  ballistic  conditions.  The  specifications  for  this  case 
are  given  by  the  first  column  of  Table  1.  The  analytic  solution  of  the  one-dimensional,  inviscid 
gas  equations  governing  the  flow  within  this  150-mm  tube  at  certain  positions  from  the  breech 
and  for  specified  times  was  obtained  by  Love  and  Pidduck.  (See  Reference  [14].)  Their  solution 
is  valid  under  the  assumptions  of  isentropic  expansion  of  each  element  of  gas,  of  constant  covo¬ 
lume,  of  an  even  integer  value  of  the  ratio  (7  +  1)/(t  -  1)  ,  where  7  is  the  ratio  of  specific 
heats,  and  the  frictionless  motion  of  the  projectile.  The  solution  is  in  terms  of  truncated  power 
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series,  and  becomes  more  complicated  as  the  number  of  rarefaction  wave  reflections  increase. 
For  the  one-dimensional  DELTA  calculation,  the  computational  mesh  (covering  the  enclosed 
cavity  behind  the  projectile)  consisted  of  four  equidistant  mesh  lines  parallel  to  the  axis  of  sym¬ 
metry  and  89  uniformly  spaced  grid  lines  orthogonal  to  the  axis  of  symmetry.  To  maintain  a 
one-dimensional  simulation  for  comparison  to  the  Love  Pidduck  results,  the  following  condi¬ 
tions 

dw^ds_=d3_  =  Q 
dr  dr  dr 


are  imposed  along  both  the  tube  wall  and  axis  of  symmetry.  The  boundary  conditions  at  the 
breech  and  projectile  are  no-slip  velocity  and  adiabatic  walls.  Love  and  Pidduck  developed  their 
solution  with  the  assumption  of  an  inviscid  isentropic  flow  which  allowed  them  to  use  a  special 
form  of  the  Noble-Abel  equation  of  state,  namely 


=  constant, 


where  subscript  zero  indicate  their  initial  values.  However,  in  the  DELTA  simulation  viscous 
effects  are  included  and  the  general  form  of  the  Noble-Abel  equation  of  state.  However,  the 
special  form  of  the  equation  of  state  used  by  Love  and  Pidduck  is  maintained  to  within  less 
than  two  percent  in  the  DELTA  simulations.  Thus,  the  non-isentropic  and  viscid  effects 
included  in  the  general  framework  of  DELTA  are  minor  for  this  one-dimensional  flow,  and  the 
comparison  of  the  analytic  solution  and  numerical  results  is  reasonable. 

The  comparison  of  the  pressure  histories  at  the  projectile  base  is  given  in  Figure  1.  The 
change  in  slope  in  the  pressure  curve  is  due  to  the  first  reflection  of  the  rarefaction  wave  at  the 
projectile  base.  The  magnitude  of  the  slope  discontinuity  of  the  pressure  curve  decreases  with 
time  due  to  the  equilibration  of  the  pressures  during  the  gas  expansion.  Another  slope  change 
exists  theoretically  at  7.137/m,  although  it  cannot  be  detected  even  in  the  graph  of  the  analyti¬ 
cal  results.  Figure  2  is  similar  to  Figure  1  except  that  the  pressures  at  the  breech  are  com¬ 
pared.  The  effect  of  the  first  arrival  of  the  wave  at  the  breech  is  more  obvious  in  both  solu¬ 
tions. 

Figures  3  and  4  show  comparisons  of  the  histories  of  the  projectile  velocity  and  projectile 
displacement  from  the  breech,  respectively.  The  large  values  of  the  tangent  to  the  curve  in  fig¬ 
ure  3  indicate  the  extreme  acceleration  the  projectile  experiences.  The  agreement  of  the  results 
show  that  the  numerical  solution  of  the  partial  differential  equations  governing  the  gas  motion 
using  the  DELTA  algorithm  is  correctly  coupled  to  the  proper  solution  of  the  projectile  motion 
Comparisons  of  the  pressure  profiles  from  the  breech  to  the  projectile  at  specified  times  arc- 
given  in  Figures  5  and  6.  The  ranges  of  pressure  values  on  the  ordinate  are  the  smallest  possi¬ 
ble  to  provide  accurate  comparisons.  Figure  5  shows  a  comparison  of  the  pressure  values  at 
2.898/m,  that  is,  after  the  rarefaction  wave  has  been  reflected  from  the  breech,  then  the 
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projectile  base,  and  is  approximately  halfway  between  them.  The  DELTA  calculation  differs  by 
0.6%  at  most  from  the  analytic  solution  values.  The  slope  discontinuity  is  smeared  out  in  the 
numerical  calculation.  Figure  6  shows  the  pressure  profiles  at  10.23  mr  which  is  near  tube-exit 
time  of  the  projectile.  The  maximum  discrepancy  between  the  two  results  is  approximately 
1.2%.  Because  the  analytic  solution  is  a  truncated  power  series  solution,  it  is  difficult  to  deter¬ 
mine  which  values  are  more  accurate. 

Next  we  compare  a  DELTA  simulation  with  time-resolved  measurements  of  the  axial  velo¬ 
city  field  inside  a  tube  behind  a  slowly  accelerating  projectile  obtained  by  laser-Doppler 
anemometry.  The  experiment  was  performed  at  Imperial  College,  London  under  funding  by 
European  Research  Office  of  the  Army  and  The  Ballistic  Research  Laboratory,  and  is  reported 
in  Reference  [15].  This  experiment  is  important  to  the  development  of  DELTA  because  it  pro¬ 
vides  the  first  transient,  two-dimensional  measurements  of  a  quantity  to  which  the  results  of  a 
DELTA  simulation  can  be  compared.  The  schematic  of  the  apparatus  is  <;  ven  in  Figure  7. 
Nitrogen  is  the  gas,  and  the  other  characteristic  of  the  experiment  are  given  in  the  second 
column  of  Table  1. 

Three  important  conclusions  of  this  study  are:  The  maximum  intensity  level  of  the  tur¬ 
bulence  was  approximately  four  percent  which  implies  a  very  low  level  of  turbulence,  the  tube 
wall  boundary  layer  remained  laminar,  and  the  heat  transfer  to  the  tube  wall  was  minimal. 
Consequently,  a  laminar  flow  simulation  with  adiabatic  walls  should  approximate  this  experi¬ 
ment.  The  t  wo-dimensional  computational  grid  had  33  uniformly  distributed  points  in  the  axial 
direction,  and  19  nonuniformly  distributed  points  in  the  radial  direction.  The  radial  grid  was 
such  that,  while  19  points  spanned  the  distance  from  the  axis  of  symmetry  (r  =  0)  to  the  wall 
(r  =  38.3 5mm  ),  12  points  were  distributed  from  r  =  32 mm  to  the  wall  and  5  points  were  dis¬ 
tributed  from  r  =  37 ,85mm  to  the  wall.  The  maximum  and  minimum  distance  between  the 
grid  points  were  6.5  mm  and  70/im,  respectively.  Constant  values  were  used  for  the  coeffi¬ 
cients  of  viscosity  and  thermal  conductivity  of  nitrogen,  namely,  17. 01 n(Pa- s )  and 
0.02524 W/(mK),  respectively.  Because  the  experimental  apparatus  was  mounted  vertically,  the 
equation  for  the  projectile  motion  (2.3)  was  changed  to  include  its  acceleration  due  to  gravity. 
Because  the  total  retarding  force  (Fp  +  fjj  +  FA  in  (2.3)  )  experienced  by  the  projectile  was  not 
determined  by  the  experiment,  no  values  of  these  forces  could  be  assigned.  Thus,  a  total 
retarding  force  profile  versus  axial  displacement  was  obtained  so  that  the  axial  velocity  of  the 
projectile  determined  by  the  DELTA  code  matched  the  experimental  values  as  shown  in  Figure 
8.  Because  the  projectile  velocity  values  agree,  so  must  the  projectile  displacement  values.  Fig¬ 
ure  9  compares  the  axial  velocity  profiles  along  the  axis  of  symmetry  at  various  times.  Both  the 
DELTA  and  experimental  results  show  a  linear  profile  from  the  zero  value  at  the  breech  to  the 
value  of  the  projectile  velocity  for  each  time.  The  axial  velocity  histories  at  76.7mm  from  the 
breech  and  at  0.5,  1.0,  2.0,  and  3.0mm  from  the  tube  wall  are  compared  in  Figure  10.  The 
values  from  the  calculation  are  within  the  scatter  of  the  experimental  data  at  radial  positions  of 
1,  2  and  3 mm  from  the  wall.  However,  the  discrepancy  between  the  values  at  the  0.5mm  posi¬ 
tion  increases  with  time.  The  same  quantities  are  graphed  in  Figure  11  but  at  153.4mm  from 
the  breech.  The  comparisons  in  Figure  11  show  similar  behavior  to  that  in  Figure  10,  but  with 


the  discrepancy  at  the  0.5 mm  position  considerably  larger.  After  a  discussion  with  the  experi¬ 
mentalists,  it  was  agreed  that  these  most  difficult  measurements  at  0.5mm  from  the  tube  wall 
are  likely  to  contain  errors  and  should  be  redone.  We  are  presently  awaiting  accurate  measure¬ 
ments  in  the  sub-millimeter  range.  However,  the  present  agreement  between  the  experimental 
measures  and  calculations  are  encouraging. 

4.  SUMMARY 


The  numerical  algorithm  encoded  in  the  DELTA  computer  code,  and  comparisons  of  its 
calculations  to  an  analytic  solution  and  experimental  measures  are  described. 

A  numerical  algorithm  to  solve  the  two-dimensional,  axisymmetric,  unsteady,  finite  volume, 
weighted  averaged  two-phase  equations  which  govern  certain  flows  inside  a  gun  tube  is  dis¬ 
cussed.  These  equations  are  in  their  most  general  form.  In  particular,  when  the  flow  regime  is 
governed  only  by  convection  and  pressure  forces,  this  general  form  automatically  gives  the  solu¬ 
tion  of  the  corresponding  inviscid  equations.  When  in  the  boundary  layer  regime,  this  general 
form  gives  without  any  assumptions,  the  solution  in  the  boundary  layer  where  viscous  forces 
dominate.  Moreover,  this  approach  naturally  provides  all  the  coupling  between  different 
phenomena  because  only  one  set  of  equations,  which  govern  all  the  phenomena,  is  solved. 
Thus,  phenomena  that  is  controlled  by  basically  inviscid  flow  but  exists  because  of  viscous 
forces,  like  the  additive  particle  laden  gun  tube  wall  boundary  layer  which  governs  heat  transfer 
to  the  gun  tube,  can  be  studied  for  the  first  time  without  assumptions  on  the  natures  of  the 
core  flow  or  boundary  layer,  the  validity  of  heat-transfer  correlations,  and/or  the  intra-flow  cou¬ 
pling  between  various  regimes. 

Two  examples  of  calculations  with  the  DELTA  code  are  presented.  These  computations  are 
limited  to  one-phase  expansion  flows  in  an  adiabatic  tube  so  to  allow  comparisons  with 
independently  obtained  data.  More  realistic  calculations  that  involve  ballistic  environments, 
heat  transfer  to  the  gun  tube  wall,  and  other  phenomena  are  presented  in  Reference  [16], 
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Figure  2.  The  comparison  of  the  pressure  histories 
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The  comparison  of  the  pressure  profiles  from  the 
breech  to  the  projectile  at  time  10.23ms. 
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Figure  7.  A  schematic  of  the  Bicen-Whitelaw  experimental 
apparatus 
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Figure  8.  The  comparison  of  the  axial  velocities  and 
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Figure  9.  The  comparison  of  the  axial  velocity  profiles  along 
the  axis  of  symmetry  at  various  times: 

(a)  4.8ms  (c)  22.8ms 

(b)  11.6ms  (d)  33.8ms 
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Figure  10.  The  comparison  of  the  axial  velocity  histories  at 
77.6mm  from  the  breech. 
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Figure  11.  The  comparison  of  the  axial  velocity  histories  at 
153.4mm  from  the  breech. 
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